#!/usr/bin/python

import sys
import getopt
from unique import unique

#-- globals --
cor = None
sen = None
ppv = None
fp = None
fn = None
tp = None
tn = None
#method = 'noencode'
method = 'encode'

def usage():
    print '''Usage : mcc.py [options] [file-real] [file-predictions]

    Options:
    --------
    -h, --help                   Display this text and exit.
    -d, --debug                  Show debugging info.
    -m, --mcc                    Show Matthews Correlation Coefficient(MCC).
    -s, --sen                    Show sensitivity(SEN).
    -p, --ppv                    Show Positive Predictive Value(PPV).
    '''

def encode(seq):
    base_map = {'A': 0, 'C' : 1, 'G' : 2, 'T' : 3, 'a':0, 'c':1, 'g':2, 't':3}

    ans = 0

    length = len(seq)

    for i in xrange(length):
        if not base_map.has_key(seq[i]): sys.exit("Wrong base type.")
        ans += base_map[seq[i]]* (4 ** (length-i-1))

    return ans

def mcc(real_file, predictions_file, option):
    global cor
    global sen
    global ppv
    global fp
    global fn
    global tp
    global tn
    right = 0
    
    file1 = open(real_file)
    lines1 = file1.readlines()
    
    file2 = open(predictions_file)
    lines2 = file2.readlines()
    
    for x in xrange(len(lines1)):
        lines1[x] = (lines1[x]).strip()
        lines1[x] = (lines1[x]).lower()
        
    for x in xrange(len(lines2)):
        lines2[x] = (lines2[x]).strip()
        lines2[x] = (lines2[x]).lower()

    lines1 = unique(lines1)
    lines2 = unique(lines2)

    # encoded lines
    en_lines1 = []
    en_lines2 = []

    if (method == 'encode'):
        for line in lines1:
            en_lines1.append(encode(line))
        for line in lines2:
            en_lines2.append(encode(line))
    
    fp = 0
    fn = 0
    tn = 0
    tp = 0

    
    #tp = len(lines1)
    #tn = 4 ** len(lines1[0]) - tp

    if (method == 'encode'):
        for i in en_lines2:
            if i not in en_lines1:
                fp += 1
            else:
                right += 1
                tp += 1

        for j in en_lines1:
            if j not in en_lines2:
                fn += 1
    else:
        for i in lines2:
            if i not in lines1:
                fp += 1
            else:
                right += 1
                tp += 1

        for j in lines1:
            if j not in lines2:
                fn += 1

    tn = 4 ** len(lines1[0]) - tp - fp -fn

    mcc = ((tp*tn)-(fp*fn))/(((fp+tp)*(tp+fn)*(tn+fn)*(tn+fp))**0.5)

    cor = mcc
    
    if option in ("-s", "--sen"):
        return (tp*1.0)/(tp+fn)
    elif option in ("-p", "--ppv"):
        return (tp*1.0)/(tp+fp)
    elif option in ("a", "--fptp"):
        return (fp, right)
    else:    
        return mcc

if __name__ == '__main__':
    try:
        opts, args = getopt.getopt(sys.argv[1:], "hdmsp", ["help", "debug", "mcc", "sen", "ppv"])
    except getopt.GetoptError:
        usage()
        sys.exit(2)

    if len(args) == 2 : mcc(args[0], args[1], None)

    for opt, arg in opts:
        if opt in ("-h", "--help"):
            usage()
            sys.exit()
        elif opt in ("-d", "--debug"):
            if len(args) != 2: usage(); sys.exit(1)
            print "FP=%d, FN=%d, TP=%d, TN=%d, MCC=%f, SEN=%f, PPV=%f , SPECIFICITY=%f" % (fp, fn, tp, tn, cor, (tp*1.0)/(tp+fn), (tp*1.0)/(tp+fp), (tn*1.0)/(tn+fp))
            sys.exit(0)
        elif opt in ("-m", "--mcc"):
            if len(args) != 2: usage(); sys.exit(1)
            print cor
            sys.exit(0)
        elif opt in ("-s", "--sen"):
            if len(args) != 2: usage(); sys.exit(1)
            print (tp*1.0)/(tp+fn)
            sys.exit(0)
        elif opt in ("-p", "--ppv"):
            if len(args) != 2: usage(); sys.exit(1)
            print (tp*1.0)/(tp+fp)
            sys.exit(0)
        else:
            usage()
            sys.exit("**** Need two files as arguments.")

    if len(args) in (0, 1):
        usage()
        sys.exit(2)

    print "MCC = %f" %   cor
